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In this paper we present a fully relativistic study of axisymmetric magnetohydrodynamic Bondi- 
Hoyle accretion onto a moving Kerr black hole. The equations of general relativistic magnetohydro- 
dynamics are solved using high resolution shock capturing methods. In this treatment we consider 
the ideal MHD limit. The parameters of interest in this study are the adiabatic constant F, the 
asymptotic speed of sound cj°, and the plasma beta parameter /3p. We focus the investigation on 
the parameter regime in which the flow is supersonic, or when > c^°. In some cases, subsonic 
asymptotic flows are considered for comparison purposes. We study the accretion rates of the total 
energy and momenta, as well as the hydrodynamic energy and momentum accretion rates. The 
models presented in this study exhibit a matter density depletion in the downstream region of the 
black hole which tends to vacuum (po = 0) in convergence tests. This feature is due to the presence 
of the magnetic field, more specifically the magnetic pressure, and is not seen in previous purely 
hydrodynamic studies. 

PACS numbers: 



I. INTRODUCTION 

The Bondi-Hoyle-Lyttleton accretion process was 
originally studied in 1942 0]. This model of the accretion 
process assumes a massive point particle travels through 
a perfect fluid background (no heat transport) with a 
Newtonian treatment of gravity. Such a configuration is 
thought to approximate the dynamics of two bodies in 
a common envelope, or the dynamics of a body travel- 
ling through the medium accreting onto an active galac- 
tic nuclei [5]. In such astrophysical systems we are able 
to safely neglect viscosity in the fluid treatment due to 
the length and velocity scales involved in the dynamics. 
However, past investigations of this phenomenon rarely 
account for the presence of magnetic fields. In many cases 
this is due to the relative simplicity of the hydrodynamic 
models. Astrophysically, it is expected that material ac- 
creting onto, say, an AGN, will be highly ionized, and 
consequently a magnetic field will be generated by the 
accreting material [7J. Due to the high conductivity of 
the accretion disk, the magnetic field is bound to the 
disk, by an effect known as flux freezing [23 . The flux 
freezing is a characteristic of ideal magnetohydrodynam- 
ics. The typical magnetic field strengths for AGN have a 
wide range range, from 10 4 G [H] to 10 10 G [22] (as cited 
by Tyul'bashev [37J). The magnetic fields observed in 
common envelopes are on the order of 10 7 — 10 8 G [33]. 

We investigate a modified version of this problem 
where the massive point particle is replaced with a black 
hole with a non-trivial radius. Consequently, we must 
consider a relativistic treatment of gravity. This line of 
research began with Petrich et al. [3T] where they found 
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a closed form solution of an ultrarelativistic spherically 
symmetric black hole and found that the flows will be 
steady. Further work by the same authors, [30], solved 
the fluid equations of motion using numerical methods. 
In [SU] they studied both the point mass model with New- 
tonian gravity, as well as a model in which the point 
mass is replaced with a spherically symmetric black hole 
and relativistic treatment of gravity. Both the Newto- 
nian and relativistic models were studied using axisym- 
metry. In both works they found that the flows were 
steady to long term evolution. Font et al. [13 EH [H] ex- 
tended this research to include an axisymmetric black 
hole, and in [TTJ they modelled a infinitely thin disk 
model, where accretion is assumed to occur only in the 
equatorial plane. In all previous research only a hydrody- 
namic background fluid was studied, and all relativistic 
flows were determined to reach a steady state. This prob- 
lem continues to be of interest in modern research such as 
in the work by Farris et al. [11] who found an application 
of the three dimensional relativistic Bondi-Hoyle accre- 
tion in binary neutron star mergers, in perturbations of 
fully three dimensional hydrodynamic models for Bondi- 
Hoyle-Lyttleton accretion [3J, and most recently in the 
search for QPO's once again using the infinitely thin disk 
approximation [35] , 

Our contribution to this research determines the phe- 
nomenology of the same system when including a back- 
ground magnetic field. We extend the original analysis 
by Font et al. [TSJ [TBI I2S] by introducing an asymptoti- 
cally uniform magnetic field as described by Wald [35]. 
In this paper, we use the ideal magnetohydrodynamic 
stress-energy tensor which reflects the presence of an em- 
bedded magnetic field. To parameterize the magnetic 
field, we use a method from plasma physics, which intro- 
duces the plasma beta parameter, (3p [T8J [T9] [25]. The 
plasma beta parameter is a ratio of the magnetic pressure 
to the thermodynamic pressure. The ideal MHD approx- 
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imation is only valid for systems in which the plasma 
beta parameter obeys the condition /3p < 1. We use 
initial parameters for our system such that the plasma 
beta parameter is less than or approximately one over 
the entire domain of integration. As reported in |16j . 
the interesting dynamics are in what is referred to as the 
supersonic regime, where the asymptotic velocity of 
the central body is greater than the asymptotic speed of 
sound cf in the fluid. 

During this study we restrict our attention to the "hot" 
relativistic equation of state, P = (T — l)pae. In partic- 
ular, we focus on two different values of the adiabatic 
constant T, one in the non-relativistic regime T = 5/3, 
and the other the relativistic equation of state T = 4/3. 
Although other values in between these are possible and 
considered physically valid, we feel that the interesting 
dynamics for a first approach to solving the GRMHD 
Bondi-Hoyle problem are captured by studying the ex- 
tremes. In a future study we will study a wider range of 
the parameters. 

One of the difficulties in studying magnetized fluids in- 
volves the enforcement of the V • B = constraint. The 
numerical treatment of magnetic fields has only become 
tractable with the techniques developed such as the con- 
straint transport method |10j , the flux transport method 
[36] . or the hyperbolic-divergence cleaning method [8]. 
In this study we implement the hyperbolic divergence 
cleaning method as described by Palenzuela et al. [2"8] . 

The outline of the paper is as follows: section [TT] will 
discuss the coordinates used in the problem. Section [TTT] 
will discuss the equations of motion used for the rela- 
tivistic ideal magnetohydrodynamic system. Section IV 
will discuss the initial conditions and boundary condi- 
tions used to perform the simulations. Section [V] will 
briefly cover the numerical methods developed and used 
for the simulations. Section I Villi will discuss the flow 
morphology and in section |VT| we draw our conclusions. 

In the rest of this paper we use geometric units, where 
G = c = I with c the speed of light in vacuum, and G 
being Newton's gravitational constant. Further we will 
use the notation that Greek scripts run over the entire 
spacetime, i.e. /i = 0,1,2,3, and Latin scripts run over 
the spatial coordinates, i.e. i = 1,2,3. All variables are 
assumed to be functions of both the radial and polar 
coordinates. 



II. COORDINATES 

In our study we are interested in the flow around a ro- 
tating black hole. Thus the line element defined by our 
black hole will be described using the Kerr spacetime. 
The original formulation of this line element contains co- 
ordinate singularities as is described in for example [6]. 
To circumvent the coordinate singularity, we use what are 
know as Kerr-Schild coordinates, as described by Font et 
al. [33]. In the Kerr-Schild coordinate system, the rotat- 



ing black hole line element is written 
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(2) 



a is the dimensionless measure of the rotation rate of the 
black hole and is related to the angular momentum of 
the black hole via J = M 2 a, where M is the mass of the 
black hole. For the present study we set M = 1 without 
loss of generality. 

In this work, we consider an axisymmetric spacetime 
geometry. As a result of this symmetry we restrict our 
study to an asymptotically uniform magnetic field which 
is aligned with the axis of rotation of the axisymmetric 
black hole. This particular magnetic field configuration 
is the only configuration that is compatible with the sym- 
metries imposed. For more general magnetic field config- 
urations, such as those presented in Bicak et al. [5], we 
will require the use of a three dimensional code, currently 
under investigation. 



III. EQUATIONS OF MOTION 

To obtain the equations of motion for the ideal rela- 
tivistic magnetohydrodynamic (MHD) system, we con- 
sider the conservation of the baryon density, the stress 
energy tensor and the Maxwell equations as seen in [§]. 
To close the system of equations, we use an equation of 
state to relate the internal energy density to the thermo- 
dynamic pressure. To decompose our system of equations 
into the 3+1 formalism we use the ADM 3 + 1 variables 
with the spacetime metric 



ds 2 



-(a 2 - /3 l f3i)dt 2 + 2f3 l dtdx i + -, ,<!.,■■ ,l.r' (3) 



where a is the lapse function, f3 l is the shift vector, and 
jij is the induced metric on the spacelike hypersurfaces. 



A. Hyperbolic Divergence Cleaning 

The evolution of the magnetic field has a physical con- 
straint V • B = 0, we have three evolution equations and 
one constraint equation. This leaves us with an over de- 
termined set of equations. Traditionally, the numerical 
treatment of a constrained system uses free evolution, 
where the evolution equations are used to evolve the sys- 
tem of equations and there is an implicit assumption that 
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the constraint will be maintained. However, when using 
free evolution, any numerical errors that arise are often 
linked to constraint violations, thus we need a method 
that enforces the constraint as the flow evolves. 

To maintain the magnetic field constraint we use the 
hyperbolic divergence cleaning method as originally pro- 
posed by Dedner and used by other groups such as 
Neilsen et al. [TJ [5] , and Palenzuela et al. [55J • 

To implement the diffusive hyperbolic method, we add 
a divergence term acting on an auxiliary field ip of the 
form V^g^tp) to the Maxwell equations V p *F ftl/ = 0, 

V„(*F^+g^)=0 (4) 

We also add a diffusive term, nn^ij) to Q to damp out 
any V • B = violations [55J, where k is a tunable pa- 
rameter. Our final expression for the Maxwell equations 
becomes 



V Jtl ( *F^ V + g^ijj) = Kn^ijj 



(5) 



In the absence of any V • B = violations we expect 
that the extra parameter ip will reduce to zero and thus 
we recover the original formulation of the ideal MHD 
equations as described in [6] for example. 



B. Stress-Energy Tensor 

We use the now standard form of the ideal magneto- 
hydrodynamic stress-energy tensor as described in, for 
example, Noble et al. [57] 

= { Po h + b 2 )u"u v + (P + b 2 /2)g» u - Wb v (6) 

where h is the enthalpy of the system 
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C. Plasma Beta Parameter 

We introduce an expression to allow us to parameterize 
the magnetic field relative to the hydrodynamic contri- 
butions. The relativistic definition of the plasma beta 
parameter [T5J[TS] is 

9P 

$p = |r- (io) 

When the magnetic field strength increases this param- 
eter decreases, and when the magnetic field strength de- 
creases this parameter increases. The magnetic fields in 
our simulations will be initialized using the asymptotic 
plasma beta parameter, /3|P, which will fully specify the 
magnetic field strength. 



D. Equations of Motion 

We use the Valencia formulation for ideal MHD. These 
are most readily found in [TH [21] thus we will not restate 
them here. We introduce 4 new variables relating to the 
divergence cleaning we implement, which are explained 
below. 



New Variables 

Considering our hyperbolic divergence cleaning 
method we have two additional conserved quantities, the 
conserved divergence cleaned magnetic density, Tl k , and 
the divergence violation field, 4", 

II' Zv' • >'«l< (11) 

* = ip/a. (12) 

We also have the corresponding primitive variables, 

B' J = W - ^'tf (13) 

■>jj = a^. (14) 

Just as in the standard GRMHD model, there are no 
known closed form solutions to the inverse relations, so 
we use numerical methods to perform this conversion 
when necessary. Primitive variable recovery is performed 
using a modified version of Mignone's one parameter in- 
version scheme [24 . This was chosen due to its simplicity, 
and when we compared this inversion scheme against the 
two variable solver promoted by Noble et al. [57] there 
was no performance difference. The scheme is readily 
found in [53] and is not repeated here. The primitive re- 
covery for variables B J and ip are trivial and are seen in 
Eqn. ([Til). 



Modified Equations of Motion 

To determine the equations of motion for the new vari- 
ables we take the 3 + 1 projection of Eqn. ([5]), to get 

dt^W + ^-yf^g{BH l -B l v 3 + ag 1 ^) 




a a 



(15) 



(16) 



The remaining equations of motion are as presented in 

El. 



INITIALIZATION AND BOUNDARY 
CONDITIONS 



IV. 



We use the method described by Font et al. [HI [16] to 
initialize the hydrodynamic variables. Since the charac- 
teristics of the system are independent of the choice of 
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initial pressure, we use — 1 and note that this choice 
is entirely arbitrary. The velocity fields are initialized 
using the following form |16j : 



1 



-_Vno COS ( 
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= 0. 



zv^ sinf 



(17) 



One may easily verify that v 2 = v 2 

The magnetic field is initialized using Wald's solution 



y 2 



Fp V — B, 



(18) 



%ii and are the azimuthal and temporal Killing 
vectors respectively, the square brackets denote anti- 
symmctrization, and a is the spin parameter of the black 
hole. Bq is a scaling factor that determines the magni- 
tude of the magnetic field. This reduces to the initial 
conditions for the magnetic field components; 
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This configuration is divergence free, and is uniform as 
r — > oo. Due to the symmetry of our configuration we 
align the magnetic field with the rotation axis of the black 
hole. 

The domain of integration for this study is defined by 
^min < r < ''max and < 9 < 7r. r m - ln is determined in 
such a way that it will always fall inside the event horizon. 
r max was set to be sufficiently far from the event horizon 
that it would be effectively considered infinity. For some 
simulations, if the radial domain is not set to be large 
enough, unphysical waves travel back towards the black 
hole, ultimately destabilizing the system. 

The boundary conditions for this problem are broken 
into three regions, two for the radial coordinate, and one 
for the angular. 

(i) In the radial direction, near the event horizon, the 
flow is strictly absorbing, also known as outflow. In 
our implementation all scalar quantities are copied 
to the ghost cell regions and all vector quantities in 
the radial direction are linearly extrapolated to the 
ghost cells. This mixing of methods prevented un- 
physical backward traveling waves exiting the event 
horizon into the rest of the domain from destroying 
the setup. 

(ii) For the boundary at the outer radial domain, we 
use prescribed boundaries "upstream" of the black 



hole, and outflow boundaries "downstream" of the 
black hole. Upstream was determined to be any 
region with angular coordinate in the range tt/2 < 
9 < 7r, while the remaining domain is considered 
downstream. The inflow condition is a prescribed 
boundary condition where we use the same asymp- 
totic values as defined in the initialization routine. 

(iii) In the angular direction, due to the symmetries 
considered in this study, reflective boundary con- 
ditions are used, which copy all scalar variables, 
and all non-0-component variables to ghost cells. 
This procedure inverts the sign of the ^-component 
of all vector quantities. 



V. NUMERICAL METHODS 

The equations of motion for our system take on the 
general form, 



0_ 

Of 
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dx l 



-9 S(Q), 



(22) 



where Q are the conservative variables, F* denotes the 
flux, and S are the geometric source terms. To solve this 
system, we developed our own high resolution shock cap- 
turing MHD code based on algorithms presented in the 
works [T|)J [201 [2SH2H]- These are all based on an inte- 



gral solution of Eqn. (22) which allows for discontinuities 
in each fluid variable. In the presence of a discontinuity 
these methods solve a one dimensional Riemann problem. 

We use the second order variable reconstruction MC 
limiter, with the HLL (Harten, Lax, and van Leer) ap- 
proximate Riemann solver. The system was time evolved 
using second order Runge-Kutta integration. It is noted 
that, like all codes using approximate Riemann solvers, 
near extremum and discontinuities the order of accuracy 
reduces to first order. For the results presented we used 
regular grid spacing with a 400 x 160 numerical grid. 

Parallelization was performed using the PAMR infras- 
tructure developed by Frans Pretorius [32] as well as di- 
rect use of the Message Passing Interface (MPI). Sim- 
ulations were performed using the woodhen cluster at 
Princeton, USA, and the WestGrid cluster in Canada. 



VI. ACCRETION PROFILES 

In this section we describe the new and major features 
of each flow studied. This is done by presenting cross 
sections of the pressure accretion profiles. 

In Fig. |3j we see the cross section of the final state 
for model M3. The upstream profile at the boundary 
smoothly agrees with the the profile, indicating that 
boundary effects will not be an issue in these simulations. 
The effects due to the presence of a magnetic field are 
minimal in the upstream region, and the magnetic pres- 
sure is uniform except in the region closest to the black 
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TABLE I: Table of parameters used for the axisymmetric 
systems studied in this paper. We perform a small sam- 
pling of the available parameter space. Only supersonic flows, 
^oc > c^°, are investigated. The "H" models are purely hy- 
dro dynamic models. 



hole. The downstream profile reveals a different picture. 
As the black hole accretes, the magnetic field builds up, 
the effects are largest near the black hole. However, the 
effects of the magnetic field are noticeable well outside of 
the hydrodynamic accretion radius, r acc , 



M 



+ (C~)2- 



(23) 



is the speed of sound in the fluid as measured by an 
asymptotic observer, M is the mass of the accretor, and 
Vrx, is the velocity of the accretor as measured by an ob- 
server at asymptotic infinity. This radius is the relativis- 
tic extension to the Newtonian Bondi radius as explained 
in detail by Petrich et al. [3D] . Matter outside this radius 
is subsonic, whereas matter within this radius has a su- 
personic flow and will inevitably fall into the black hole. 
We refer the reader to [9] for details of the Newtonian 
calculations and more history. For the relativistic flows, 
this radius is fairly small, typically 0.5 < r acc < 1.5. For 
the sake of consistency, we choose to measure the accre- 
tion rate at the event horizon, where we can focus on 
the material as it falls into the black hole. The radial 
location of the event horizon only depends on the mass 
and the angular momentum of the black hole, so this is 
constant for each choice of fluid parameters. 

The profiles presented in Figs. [T] [2] and [3] have similar 
general appearance. The upstream pressure profile takes 
on the appearance of a stationary Bondi accretor, as was 
originally presented in Font et al. [TB] - We find that the 
similarity continues when comparing the upstream pro- 
file of our models to the magnetic Bondi accretor. Just 
as in [16] the downstream region is significantly different 
from the stationary accretor. In all cases the downstream 
profile indicates that the magnetic pressure builds up sig- 
nificantly in this region, and dominates the total pressure 
close to the event horizon of the black hole. An interest- 
ing trend occurs when modeling Ml and M3, the up- 
stream profiles have the exact same trends; however, for 



model M3 in the downstream region the magnetic pres- 
sure is more than three times larger, and the effect has a 
shorter range. The profile for model Ml shows that the 
magnetic pressure extends to approximately r = 20M , 
while the same profile for model M3 shows that the mag- 
netic pressure only extends to approximately r — 15M . 

In Fig. [4] we present a pseudocolour plot for the ther- 
mal pressure for model Ml with a contour plot of the 
magnetic pressure. The magnetic field wraps around the 
black hole, as is seen by the magnetic pressure contour 
in the upstream region that traces the edge of the Mach 
cone in the downstream region. This forces the matter to 
decrease there. The magnetic reconnection in the down- 
stream region is sufficiently weak that the magnetic field 
lines pile up in the downstream region. The magnetic 
pressure forces the downstream matter away from the 
axis of symmetry creating an evacuated region, which is 
similar to the effect as seen in the plasma depletion region 
associated with the earth's magnetosphere. 

In Fig. [5] we see the pseudocolour plot of the thermal 
pressure, with a contour plot of the magnetic pressure 
for model M2. In addition we present the magnetic field 
vectors. Outside of the Mach cone the magnetic field 
points in the same direction as the asymptotic magnetic 
field, however the field reduces magnitude, and reverses 
direction when it crosses into the tail shock. The magni- 
tude of the magnetic field is represented by the size of the 
head of the vector. The magnitude is the largest closest 
to the black hole. 

Figures [6] [7J and [8] present evidence that there is a 
balance between the thermal pressure and the magnetic 
pressure, in particular the total pressure has a similar 
profile to the hydrodynamic models. 

A general feature of all relativistic Bondi-Hoyle accre- 
tion is the presence of a Mach cone [T5HT71, [35] . This 
cone attaches to the downstream side of the black hole 
with what is known as the opening angle. We observe 
the opening angle of the Mach cone, seen in Figs. 6|7 
and [8] The opening angle is known to be a function of 
the parameters of the fluid for Newtonian systems |13j . 
and is seen in the study performed by Font et al. [To] , 

In Figs. |9j [TUJ and [TTJ we present the cross section of 
the pressure accretion profiles on r — 2M for models 
Ml, M2, and M3, as well as models HI, H2, and H3 
for comparison. The balance between the two sources 
of pressure, magnetic and thermal, is clear in these cross 
sections since the total pressure exhibits a similar form as 
those seen in purely hydrodynamic flows [16]. The effect 
of the magnetic field increases the overall pressure in this 
cross section, especially along the axis of symmetry, and 
widens the Mach cone opening angle. The presence of the 
magnetic field amplifies the downstream pressure along 
the axis of symmetry and widens the Mach cone angle 
of attachment. Since the region downstream, inside the 
Mach cone, is the location of the maximum accretion, 
we expect that the larger cross sectional area leads to an 
increased accretion rate, as is seen in Fig. [I3j In Fig. [Tl] 
we see that model H3 presents a maximum pressure on 
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FIG. 1: We show the pressure cross-section in the upstream region (left) and downstream (right) along the axis of symmetry 
plotted from Teh < r < 50 for model Ml. We see that the asymptotic behaviour of our system agrees with our boundary 
conditions, indicating that boundary effects are negligible for these simulations. We see a small effect due to the presence of the 
magnetic pressure in the upstream region; however, in the downstream region the magnetic effects are far more dominant in 
the region ?-eh < T < 10 which is outside of the expected hydrodynamic accretion radius as defined in Eqn. ( 23 1. The upstream 
thermal pressure agrees with the hydrodynamic pressure from model HI, and the downstream total pressure profile agrees 
with the pressure profile of the same model. The magnetic pressure profile was shifted by one to fit the scale presented here. 
As is described in [TH] the upstream pressure cross section takes on the appearance of a Bondi accretor while the downstream 
region is far more extreme. 




FIG. 2: We show the pressure cross-section in the upstream region (left) and downstream (right) along the axis of symmetry 
plotted from Teh < r < 50 for model M2. We see a similar profile as seen in Fig. [T] In this model the magnetic pressure 
has a much larger amplitude upstream than for model Ml, and shifting the data was not necessary. Upstream we also see 
a deviation between the hydrodynamic model and the thermal pressure. In the downstream region the magnetic effects are 
far more dominant in the region Teh < r < 20 well outside the hydrodynamic accretion radius. The downstream the total 
pressure from model M2, and the hydrodynamic pressure from model H2 are again very similar. 
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FIG. 3: We show the pressure cross-section in the upstream region (left) and downstream (right) along the axis of symmetry 
plotted from Teh < r < 50 for model M3. The profile in this model is very similar to that seen in Fig. [I] including the 
upstream similarities between the hydrodynamic pressure and the thermal pressure. The effect of the rotating black hole is 
seen by the shortened magnetic pressure profile, which decreases to teh < r < 6 which, as in Fig. [5] is outside of the expected 
hydrodynamic accretion radius. Careful inspection also shows that the downstream thermal pressure profile is dramatically 
different between models Ml and M3. The thermal pressure maximum occurs around r = 3M here whereas in model Ml the 
maximum is on the event horizon. Like in Fig. [I] the magnetic pressure was shifted by one to fit the scale. 




FIG. 4: The thermodynamic pressure (colour plot) and the magnetic pressure (contours) on a logarithmic scale for model Ml. 
We see that there is a balance between the two sources of pressure. This is made more apparent in figure[6] We also see that the 
magnetic field wraps around the black hole as shown by the dark blue contour connecting the magnetic pressure downstream 
to the magnetic pressure upstream. The magnetic reconnection downstream is sufficiently weak that the magnetic field lines 
pile up on the downstream side. 



the axis of symmetry/rotation, whereas the total pressure 
for model M3 decreases on the axis. Further investigation 
is necessary to relate the spin of the black hole to the 
pressure decrease in the magnetic field cases. 



In models Ml and M2 we see that the magnetic pres- 
sure in the centre of the Mach cone has a much larger 
amplitude than the surrounding thermal pressure, while 
in model M3 the magnetic pressure has a smaller ampli- 
tude than the surrounding thermal pressure. Clearly, by 
allowing the magnetic field to "wind" with the black hole 
the magnetic pressure decreases which reduces the open- 
ing angle of the Mach cone. Since we also see that when 
we turn on the rotation, that there is very little change 
in the opening angle between the hydrodynamic model 
and the magnetohydrodynamic model. However, we do 
need to perform a larger parameter survey to determine 
how the magnetic fields are affected by the rotation of 
the black hole. 



VII. ACCRETION RATES 

In this section, we discuss the measured quantities that 
we observe for each run. The first several are simply the 
same quantities as developed by Petrich et al. |30j and 
later used by Font et al. [TBI IT?] . The accretion rates 
come from conservation laws at asymptotic infinity. We 
calculate the mass accretion rate from the mass of the 
accreted matter, 

m = J d 3 xp (24) 

which when we take the time derivative and substitute 
the mass flux into the integral becomes, 

772 = j d 3 xd l (pv l ) = 4tt j> s /^gpv r dd. (25) 

We calculate the energy and momentum accretion rates 
in a similar fashion, 

Q (t) = -4vr f d 3 xT^n„ 1% (26) 
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FIG. 5: The thermodynamic pressure (colour plot) and the magnetic pressure (contours) on a logarithmic scale for model M2. 
We superimpose the magnetic field vectors. We see that the downstream side has a distinct region where the magnetic field 
vectors diminish and switch direction. The thermal pressure balances the magnetic pressure in this region. This is a similar 
effect as seen in the plasma depletion region associated with the earth's magnetosphere. 



where is the L-th. Killing vector of the system, and 
is a normal vector to the 3+1 hypersurface. With 
the given 3 + 1 description, we use = (—a, 0,0,0) so 
the above equation reduces to, 



)(') 



4tt / d A xT 



tv a (t) .' 



(27) 



For the axisymmetric black hole there are two Killing 
vectors, one temporal and one azimuthal. These lead to 
conservation of energy and azimuthal angular momen- 
tum. At asymptotic infinity, the metric takes the form 
of Minkowski space which is maximally symmetric |17j . 
Thus we may consider the radial Killing vector. By the 
same logic we are able to use the polar Killing vector as 



well; however, due to the symmetries involved in this 
study, no relevant physics can be extracted from this 
value. 

To obtain the momentum accretion rates, we take the 
time derivative of Eqn. (27), where the time deriva- 
tive commutes with the spatial integral. Subsequently, 
we replace the time derivative quantities with their flux 
counterpart to get an expression for the accretion rates 
[TBI 133] . Since our HRSC scheme is at most second or- 
der, we use a second order numerical integration scheme 
to solve the resulting expressions. 

The quantities we are specifically interested in present- 
ing are the radial momentum accretion rate, the energy 
accretion rate, and the mass accretion rate. 
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FIG. 6: The thermodynamic pressure (left), the magnetic pressure (middle) and the total pressure (right) for model Ml. We 
see that the total pressure takes on a familiar form, exhibiting a tail shock as seen in the work by Font et al. [16] . 




The first quantity considered is the energy accretion 
rate. The equation is stated here and is easily derived 
using Eqn. (26). 



Q 



(i) 



Q 



(*) 



E 




gaT tf d 3 x = - J ^a 2 T tt d 3 x 
jEd 3 x 
~gT tr d6 



y/^gT^T\ u drd9. (28) 

Likewise we have the momentum accretion rates 
=P r = I y/^gT rr d6 



-gT^T\ v drdO, 



(29) 



where teh is the radial location of the black hole event 
horizon and r m denotes the radius that we measure the 
accretion rate. In our measurements, we are interested in 
capturing dynamics on the event horizon, so r rn = r^n, 
thus all volume integrals above are zero. 



To normalize the mass accretion rate, previous authors 
used the mass accretion rates determined at the sonic 
point for a relativistic Bondi (stationary) accretor [16l 

M, 



M 



A-kXM 2 Poo-tub 

(<4 + (C) 2 ) 3/2 ' 



(30) 



A is a dimensionless form factor as described in Shapiro 
and Teukolsky [33 . When normalizing the radial mo- 
mentum accretion rate, they scaled the solution by both 
the mass accretion rate and the asymptotic velocity. 



ATTXM 2 p O0 m B v O 3 



(31) 



We do not use these factors for two reasons. First, we 
measure the accretion rates at the event horizon not the 
accretion radius; second, we are measuring the accretion 
rates for a magnetohydrodynamic system, and the above 
simple forms do not take into account magnetic effects. 
Instead we time evolve the relativistic Bondi-Hoyle sys- 
tem using = 0, therein solving the relativistic mag- 
netized Bondi problem, starting from a uniform density 
background. When the flow reaches a steady state, we 
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FIG. 9: A cross section at r = 2M for model Ml of the 
thermodynamic pressure, magnetic pressure, and the total 
pressure, as well as the pressure for model HI. The balance 
between the two sources of pressure is more clear in this cross 
section since the total pressure exhibits a similar form as those 
seen in purely hydrodynamic flows [16] . The effect of the mag- 
netic field increases the overall pressure in this cross section, 
especially along the axis of symmetry, and widens the Mach 
cone angle of attachment. 

use those values of the accretion rates to normalize the 
accretion rate measurements for the relativistic Bondi- 
Hoyle systems. 

Due to the symmetry of the system the azimuthal ve- 
locity field is exactly zero and, consequently, is not plot- 
ted. 

We see that in all accretion plots that the flows reach a 
steady state. We also see the general trend that all mag- 
netic flows have a wider Mach cone opening angle and 
consequently all magnetic flows experience a larger mass 
and energy accretion rate. This trend also extends to 



FIG. 10: A cross section at r = 2M for model M2 of the ther- 
modynamic pressure, magnetic pressure, and the total pres- 
sure. We also plot the same cross section of the pressure for 
model H2. The balance between the two sources of pressure 
clearer in this cross section since the total pressure exhibits a 
similar form as those seen in purely hydrodynamic flows [16] . 
As was the case in model Ml, the presence of the magnetic 
field amplifies the downstream pressure along the axis of sym- 
metry and widens the Mach cone angle of attachment. Since 
the region downstream, inside the Mach cone, is the location 
of the maximum accretion, we expect that the larger cross 
sectional area leads to an increased accretion rate, as is seen 
in Fig. [13] 



the radial momentum accretion rate for all models with 
r = 4/3; however, for models with T = 5/3 we find that 
this accretion rate decreases for larger amplitudes of the 
magnetic field. This is attributed to the relative build 
up of material in the downstream region. For all models 
studied we have an evacuated region downstream; how- 
ever, it is observed that the models with r — 4/3 have 
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FIG. 11: A cross section at r = 2M for model M3 of the 
thermodynamic pressure, magnetic pressure, and the total 
pressure, and the pressure profile for model H3. The bal- 
ance between the two sources of pressure is more clear in this 
cross section since the total pressure exhibits a similar form 
as those seen in purely hydrodynamic flows [16 . As in the 
case for model Ml, we find that the total pressure exceeds 
that for model H3, where no magnetic field is present; how- 
ever, model H3 presents a maximum pressure on the axis of 
symmetry/rotation, whereas the total pressure for model M3 
decreases on the axis. Further investigation is necessary to 
relate the spin of the black hole to the pressure decrease in 
the magnetic field cases. 



a narrower high magnetic pressure region, which leads 
to a narrower region of the shock tail that experiences 
density evacuation. As such this leaves a larger region 
of the Mach cone in contact with the black hole, and 
clearly, overall larger accretion rates. For T = 5/3 mod- 
els the magnetic pressure causes a much larger density 
depletion region. The wider magnetic pressure profile 
decreases the contact between the matter and the black 
hole, decreasing the drag between the two. The magni- 
tude of the density in contact with the black hole is larger 
than the purely hydrodynamic simulation, accounting for 
the larger mass and energy accretion rates. This is also 
interesting behaviour, since it indicates that a cold fluid 
will experience less drag than the hot fluid. 

In Figs. [12] [l3j and [14] we compare the accretion rates 
for the radial momentum, the energy and the mass. Just 
as presented in [16] , the accretion rates are dependent on 
the choice of adiabatic constant. 

Figures [15} [TH] P^l display the mass, energy and radial 
momentum accretion rates for models Ml, M4, M5, and 
HI. As we expect, as the magnetic field reduces in mag- 
nitude, the accretion rates approach the hydrodynamic 
rates. In the blown up region within each plot we see the 
accretion rate as a function of the plasma beta parameter. 



FIG. 12: The radial momentum accretion rates for models 
Ml, M2. We see that the adiabatic constant plays a significant 
role in the amount of drag the moving black hole experiences, 
which agrees entirely with previous studies. 



The energy and mass accretion rates decrease as a func- 
tion of the magnetic field strength, and converge to the 
hydrodynamic accretion rate in the zero magnetic field 
limit. Although we do note that in Fig. [15] the presence 
of even a weak magnetic field damps out the oscillatory 
accretion rate present in the hydrodynamic study. The 
oscillatory behaviour of the hydrodynamic solution is dis- 
cussed in |16[ 135] . The radial momentum accretion rates 
also increase as a function of the magnetic field strength. 

Figures [18] [T9[ [20] present the mass, energy and radial 
momentum accretion rates for models M2, M6, M7, M8, 
and H2. As we expect from previously discussed mod- 
els, as the magnetic field reduces in magnitude, the ac- 
cretion rates approach the hydrodynamic rates. In this 
case we find that the magnitude of the radial momen- 
tum accretion actually decreases with increased magnetic 
field strength. This is in contrast to the hot plasma 
case, where the drag experienced by the spherically sym- 
metric black hole increases with increased magnetic field 
strength. 

Figures [2lJ [22] [23] show the mass, energy and radial 
momentum accretion rates for models M3, M9, M10, and 
H3. These plots indicate that when accreting onto an 
axisymmetric black hole the general trends of the flow 
do not alter significantly. Just as in Figs. [15] [16] and 
[FT] the presence of a magnetic field increase all accretion 
rates. 

The simulations presented here all had the property 
that they converge to the asymptotic solution shortly af- 
ter 600M. 
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FIG. 13: The energy accretion rates for models Ml, M2. 
When comparing the Mach cone opening angles between mod- 
els Ml and M2, we see that the accretion rates are related to 
the opening angle, the larger the cross section in the shock 
cone, the more energy in accreted. This finding is in agree- 
ment with the results of the purely hydrodynamic models in 
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FIG. 14: The mass accretion rates for models Ml, M2. Just 
as in the case of drag, the mass accretion rate is significantly 
affected by the value of the adiabatic constant used. 



FIG. 15: The mass accretion rates for models Ml, M4, M5, 
and HI. We see that the mass accretion rate is only signifi- 
cantly impacted by the presence of a large magnetic field as 
simulated in model Ml. As we expect, the mass accretion 
rates decrease as a function of the magnitude of the plasma 
beta parameter. What is of particular interest is the oscil- 
latory behaviour of the purely hydrodynamic model, investi- 
gated by |35) is damped out quickly by the presence of even 
a weak magnetic field as presented in model M5. 
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VIII. FLOW MORPHOLOGY 

The results of the simulations using the parameters 
found in table [T] are discussed here. We present the fi- 
nal state of a sampling of the parameters surveyed in the 
spherically symmetric evolution. We see that the ther- 



FIG. 16: The energy accretion rates for models Ml, M4, M5, 
and HI. Just as in Fig. |15[ the energy accretion rate is only 
significantly impacted by the presence of a large magnetic 
field. Indicating that the increase in the Mach cone opening 
angle allows more fluid to accrete. 
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FIG. 17: The radial momentum accretion rates for models 
Ml, M4, M5, and HI. We see that like all the previous mea- 
surements for these models that the larger opening angle leads 
to an increased drag experienced by the star. 
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FIG. 18: The mass accretion rates for models M2, M6, M7, 
M8, and H2. We see that the mass accretion rate is only 
significantly impacted by the presence of a large magnetic 
field as simulated in model M2. As we expect, the mass ac- 
cretion rates decrease as a function of the magnitude of the 
plasma beta parameter. What is of particular interest is the 
oscillatory behaviour of the purely hydrodynamic model, in- 
vestigated by [35] is damped out quickly by the presence of 
even a weak magnetic field as presented in model M5. We 
note that for this adiabatic constant the oscillations in the 
hydrodynamic simulation do not exist. The mass accretion 
rate increases monotonically as a function of the plasma beta 
parameter. 



FIG. 19: The energy accretion rates for models M2, M6, M7, 
M8, and H2. Just as in Fig. |18| the energy accretion rate is 
only significantly impacted by the presence of a large magnetic 
field, indicating that the increase in the Mach cone opening 
angle allows more fluid to accrete. As is expected, the en- 
ergy accretion rate increases in a similar fashion as the mass 
accretion rate 1301. 



modynamic quantities establish a steady state solution. 
As in the general hydrodynamic models, we find that in 
the upstream region the contours are smooth, but in the 
downstream region there is the presence of a Mach cone. 
As with the original hydrodynamic case we also look at 
the profile along the cf> — and cf> — it lines. 

The simulations result in a rest mass evacuation im- 
mediately downstream of the black hole. Upon further 
analysis, this is the result of a buildup of the magnetic 
pressure downstream in the same region. The baryonic 
particles are transported along the magnetic field lines 
away from this region. As we see in figures [6] and [7J the 
convergence tests as seen in Figs. [26] and [27] indicate that 
the region is tending to zero rest mass density. This phe- 
nomenon is similar to the behaviour of the particles in 
the solar wind as they interact with the earth's magne- 
topause [35J ■ In that setup, [35], a region known as a 
plasma depletion layer develops in the upstream side of 
the earth moving through the solar wind from the sun. 
In said region the magnetic field strength increases and 
as a result the plasma density decreases. Although this 
process is most dominant in the upstream side of the 
earth's trajectory there is evidence that it also exists for 
the downstream side when the magnetic reconnection is 
insufficient to relieve the pile-up of magnetic field lines. 
In our system, the pile-up occurs primarily in the down- 
stream side, which is likely due to the fact that the fluid 
itself is magnetic, and that the black hole is not assumed 
to have its own magnetosphere in our model. Future 
models may include such features, but such a study would 
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FIG. 20: The radial momentum accretion rates for models 
M2, M6, M7, M8, and H2. We see that unlike all the previous 
measurements for these models that the larger opening angle 
leads to an decreased drag experienced by the star. 
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FIG. 21: The mass accretion rates for models M3, M9, M10, 
and H3. We see that the mass accretion rate is only signifi- 
cantly impacted by the presence of a large magnetic field as 
simulated in model M3. As we expect, the mass accretion 
rates decrease as a function of the magnitude of the plasma 
beta parameter. First, we see no evidence of oscillatory mass 
accretion as seen in model HI. We note that for this adiabatic 
constant the oscillations in the hydrodynamic simulation do 
not exist. The mass accretion rate increases monotonically as 
a function of the plasma beta parameter. The mass accretion 
rate quickly diverges towards the end of the simulation. This 
is determined to be a result of the resolution of the simulation. 



FIG. 22: The energy accretion rates for models M3, M9, M10, 
and H3. Just as in Fig. |21[ the energy accretion rate is only 
significantly impacted by the presence of a large magnetic 
field, indicating that the increase in the Mach cone opening 
angle allows more fluid to accrete. As is expected, the en- 
ergy accretion rate increases in a similar fashion as the mass 
accretion rate [30| . 



be purely academic, since uncharged black holes are not 
expected to have their own magnetic fields. 

Following the flow from upstream to downstream 
within a few Schwarzschild radii of the black hole reveals 
a very similar morphology as the hydrodynamic models. 
The flow is drawn to the black hole via gravitational at- 
traction. As the material flows past the black hole, it is 
attracted to the hole by gravitational forces. The angu- 
lar momentum of the fluid prevents the direct inflow of 
all the gravitationally bound material. This angular mo- 
mentum is lost as the material moves downstream and 
begins to converge on the axis of symmetry. Here the 
fluid increases in density and pressure. This is where the 
similarities end. Closer to the black hole, we see a no- 
ticeable difference. While the same flow process occurs 
the effects of the magnetic field are far more obvious. 
The downstream pressure increases, including the mag- 
netic pressure. As more fluid accumulates the frozen flux 
tubes also accumulate, which increases the local magnetic 
field strength. Eventually the magnetic pressure begins 
to dominate the local flow, and the density in that region 
begins to decrease as the fluid pressure outside of this re- 
gion deflects new fluid from entering, but at the same 
time the magnetic pressure decelerates the fluid and di- 
rects it into the black hole. A balance between magnetic 
pressure and thermodynamic pressure is established and 
the flow reaches a steady state. 

In Fig. [9] we have a cross section of model Ml. We 
see that the total pressure of the system appears contin- 
uous; however, when we look at the contributions from 
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FIG. 23: The radial momentum accretion rates for models 
M3, M9, M10, and H3. We see that unlike all the previous 
measurements for these models that the larger opening angle 
leads to an decreased drag experienced by the star. 

the different pressures the story is very different. In the 
evacuated region, we see that the magnetic pressure is 
dominant, showing that the force due to the magnetic 
field is strongest in this region. The radial velocity of 
the fluid in that region becomes negative indicating flow 
into the black hole, which is being forced faster than a 
hydrodynamic model due to the pressure from the mag- 
netic field. As the oncoming flow reaches the tail shock, 
it first interacts with a region dominated by the thermo- 
dynamic pressure. If the fluid velocity is great enough 
the fluid will pass the first shock point, along the way 
interacting with a relatively minor build-up of magnetic 
field within the thermal pressure wall. Once past the 
wall, the fluid reaches a region where the thermal pres- 
sure suddenly drops and the magnetic pressure dramati- 
cally increases. It is in this region that the magnetic field 
is at its strongest, which forces the fluid to follow the path 
of least resistance, either the fluid accretes onto the black 
hole or flowing faster upstream, either way causing a den- 
sity depletion. The low point in the magnetic pressure 
denotes the location where the magnetic field changes di- 
rection and points towards the black hole. This is seen 
most easily in Fig. [5j 

When investigating axisymmetric Bondi-Hoyle accre- 
tion onto a rotating black hole we found that the gen- 
eral flow morphology was substantially different to the 
non-rotating black hole with the same equation of state. 
Which is most dramatic when observing the plasma beta 
parameter as a function of the computational domain in 
Fig. [24J 

Finally, in Fig. [24] we present the plasma beta pa- 
rameter through the domain, after each simulation has 
reached a steady state. We see that a black hole trav- 



elling well above supersonic speeds produces a readily 
apparent tail shock. Also note the peaks that form along 
the inside edge of the tail shock for model M3. This is 
a feature only seen in simulations with a rotating black 
hole. The plasma beta parameter in the tail clearly fol- 
lows a similar form as the tail shock near the black hole; 
however, for the spherically symmetric black holes, as we 
move further away downstream from the black hole the 
plasma beta parameter begins to smoothly reach asymp- 
totic values. The rotation of the black hole clearly winds 
the magnetic field which may be seen in the maximum of 
ftp in model M3 being an order of magnitude larger than 
that seen in model Ml. 



IX. CONCLUSIONS 

In this paper we presented a preliminary survey of 
the magnetized Bondi-Hoyle accretion onto both spheri- 
cally symmetric and rapidly rotating axisymmetric black 
holes with the rotation axis aligned with the asymp- 
totic uniform magnetic field. We used a high-resolution 
shock-capturing scheme with an approximate Riemann 
solver. To maintain the V • B = constraint, we used 
a hyperbolic-divergence cleaning method and used the 
auxiliary field as more than just a means to remove di- 
vergence violations, but also as a well behaved function 
that may be used to determine the convergence rates for 
magnetized systems. 

We have established that when magnetic fields are in- 
troduced to the typical hydrodynamic calculations, the 
evolution of an axisymmetric system exhibits behaviour 
that is qualitatively broadly similar to the purely hydro- 
dynamic counterpart in [15j [16] but differs quantitatively 
in a range of crucial ways. Ultimately, we have shown 
that all parameters surveyed resulted in a steady state so- 
lution with an evacuated region downstream of the black 
hole. This evacuation is a result of a build-up of mag- 
netic pressure in the same region, which forces all matter 
out of this region. Future studies will investigate a wider 
range of parameters and will attempt to match the pa- 
rameters investigated to astrophysical observations. We 
have also seen evidence that although a rotating black 
hole in a hydrodynamic simulation aligned in a config- 
uration similar to ours yields very little in the way of 
different morphology, the introduction of an ideal mag- 
netic field is sufficiently different to make the dynamics 
much more interesting. 

The presence of the magnetic field is seen to increase 
the Mach cone opening angle, and consequently increase 
the accretion rates for both the mass and energy. How- 
ever, we see that the radial momentum accretion rates, 
or drag, decreases as a function of the magnetic field 
strength for cold plasmas, while hot plasmas experience 
more drag as a function of the magnetic field strength. 
This is attributed to the fact that the magnetic pressure 
of a r = 5/3 fluid downstream of the black hole is wider 
(relative to the T = 4/3 models), and creates a wider 
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FIG. 24: Plasma beta parameter, /3p(r,8) for the axisymmetric evolution for model Ml (left, t = 2500M), M2 (centre, 
t — 1500 M), and M3 (right, t = 2500M). These snapshots were taken well after the flow reaches its steady state. We see that 
a black hole travelling well above supersonic speeds produces a readily apparent tail shock. Also note the peaks that form 
along the inside edge of the tail shock for model M3. This is a feature only seen in simulations with a rotating black hole. 
The plasma beta parameter in the tail clearly follows a similar form as the tail shock near the black hole; however, for the 
spherically symmetric black holes, as we move further away downstream from the black hole the plasma beta parameter begins 
to smoothly reach asymptotic values. The rotation of the black hole clearly winds the magnetic field which may be seen in the 
maximum of f3p in model M3 being an order of magnitude larger than that seen in model Ml. 



evacuated region, and less surface area of the black hole 
event horizon is exposed to the flow in a region where all 
accretion is expected to be maximal. 

This was a phenomenological study, a detailed analysis 
will require the use of less diffusive Riemann approxima- 
tions such as the Roe or Marquina solvers as used in re- 
cent hydrodynamic studies of the Bondi-Hoyle accretion 
problem in [35]. Further, a detailed study of the param- 
eters needed in the hyperbolic diffusion cleaning method 
will need to be refined. Future studies will also involve 
a full three dimensional code which will allow a general 
asymptotic magnetic field configuration as prescribed by 
Bicak [5J, in order to capture more astrophysically justi- 
fiable configurations. This line of study is far from com- 
plete, as this configuration will provide an excellent and 
simple testbed for more general fluid models as the field 
grows increasingly more complicated. 
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X. APPENDIX: CODE VERIFICATION 

To verify the numerical accuracy of the code we per- 
formed convergence testing of two different quantities. 
First we used a global quantity, the L2-norm of the auxil- 
iary variable 'F, and second we used a cross section of the 
pressure accretion profile. We considered a cross section 
of the data at a constant radial coordinate, and indepen- 
dently a slice through a constant value of the angular- 
coordinate. Different techniques for convergence testing 
HRSC schemes is an active research avenue and there is 
no universally accepted scheme at the time of writing. 
Convergence tests for models Ml and M3 may be seen 
in Figs. [25] [26] and [27} While we see that the deple- 
tion region has not converged by our selected resolution, 
the various accretion rates are within convergence. Fu- 
ture work will investigate the different parameters using 
higher resolution to work in the convergent regime of the 
depletion region. 

To determine the convergence of the divergence clean- 
ing variable ip to zero, we calculate the L2-norm of this 
scalar field. While this would at best indicate a first order 
convergence rate, the fact that this term converges to zero 
is critical in showing the validity of our selected method 
for handling divergence violations. A sample convergence 



test for the tp variable is seen in Fig. 25 



The disagreement between the true second order con- 
vergence value and our system comes from the special 
treatment of shock regions and local maxima where the 
system reduces to first order convergence as mentioned 
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FIG. 25: The convergence of the L2 norm of the auxiliary 
field ip for model M3. It is clear that the system is convergent. 
Level 1 uses a grid resolution 200 x 80, level 2 denotes 400 x 
160, and level 3 denotes 800 x 320. 
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FIG. 26: The convergence of the rest mass density density 
po cross section at r = 2M for model Ml. We see that the 
evacuated region near the axis of symmetry is converging to 
zero. Level 1 uses a grid resolution 200 x 80, level 2 denotes 
400 x 160, and level 3 denotes 800 x 320. 
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evacuated region near the axis of symmetry is converging to 
zero. Level 1 uses a grid resolution 200 x 80, level 2 denotes 
400 x 160, and level 3 denotes 800 x 320. 
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